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Abstract. Astronomers have constructed models of globular clusters for over 100 years. These 
models mainly fall into two categories: (i) static models, such as King's model and its variants, 
and (ii) evolutionary models. Most attention has been given to static models, which are used 
to estimate mass-to-light ratios and mass segregation, and to combine data from proper mo- 
tions and radial velocities. Evolutionary models have been developed for a few objects using 
the gaseous model, the Fokker-Planck model, Monte Carlo models and Af-body models. These 
models have had a significant role in the search for massive black holes in globular clusters, for 
example. 

In this presentation the problems associated with these various techniques will be summarised, 
and then we shall describe new work with Giersz's Monte Carlo code, which has been enhanced 
recently to include the stellar evolution of single and binary stars. We describe in particular 
recent attempts to model the nearby globular cluster M4, including predictions on the spatial 
distribution of binary stars and their semi-major axis distribution, to illustrate the effects of 
about 12 Gyr of dynamical evolution. We also discuss work on an approximate way of predicting 
the "initial" conditions for such modelling. 

Keywords, methods: numerical, globular clusters: general, globular clusters: individual (M4) 



1. Introduction 

1.1. Some astrophysical questions 

There are many reasons for constructing a dynamical model of a globular cluster, but 
they fall into two broad categories. First there are problems that can be tackled by 
constructing static (equilibrium) models, such as 

(a) Inferring the mass from the surface brightness profile, radial velocities, proper 
motions and mass functions. In this way one can estim ate the total mass o f stars below 
the observational limit, such as faint white dwarfs fe.g. iDrukier et al.lll988l ) 



(b) Inferring the global mass function from local mass functions (e.g. iRicher et al 
120041 ): then one can address the question of whether this is the same for all clusters. 

(c) Measuring cluster distances by comparison of radial velocities and proper motions: 
the model is used to correct for rotation, or to link the differ ent locations and stella r 
components which are observed by the different techniques (e.g. Ivan de Ven et al.|[2006h . 

Then there is a second range of questions which require dynamic evolutionary models, 
i.e. questions such as 

(a) Inferring the primordial mass function from local, present-da y mass functions: the 
mode l is used to correct for preferential escape of low-mass stars (e.g. lBaumgardt fc Makino 

ill). 

(b) Inferring primordial parameters of the binaries from the ir present-day stati stical 
properties: primordial abundance, period distribution, etc (e.g. iKroupa et al.ll200ir ) 
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(c) Determining the effect of dynamics on th e estimation of m ass through the virial 
theorem, which is affected by mass segregation ( Fleck etld1l2006l) 

Note that these last few references do not deal with specific objects (which is the focus 
of the rest of this review) but with general trends. 

This paper begins with a review of the methods and observational constraints which 
have been used to construct models of individual globular clusters, often with a view 
to answering the above types of question. Then we focus on one particular method, the 
Monte Carlo method, and its application to the nearby globular cluster M4. 

1.2. Methods and constraints 

1.2.1. The methods 

A number of techniques have been used to construct static models, to answer questions 
of the first type: 



Plummer's m odel dPlummer 1911 ) 

King's model ( King 19661: Peterson fc Kind 



975) 



1963) 



(a) 
(b) 

(c) Anisotropic models (Michic & Bodcnhcimcr _ 

(d) Multi-mass models dGunn fc Griffiij|l979t IPryor et al.lll986t iDubath et al.|[l990h 

(e) Non-parametric models (Gebhardt & Fischer 1995) 
(/) Schwarzschild's m ethod 
(g) Jeans' equations ([Leonard et al. 



van de Ven et al 
1992fi~ 



2009) 



Even this list may not be exhaustive. 

Much less work has been done on dynamical evolutionary models of specific globular 
clusters, but the method s include 

(a) Gas/fluid models (jAngeletti et al.lll980l [M3]) 

[&) Fokke r-Planck models (Co lin and co-workers: lGrabhorn et al.lll992l [N6624] jDull et al 

119971 [M15]: lDrukienll993l Il995l [N6397]: Phinney 1993 [M15]) 

(c) Monte Carlo models (Giersz fc Heggie 2003 [u> Cen]) 

(d) TV-body models 

The last of these should really not be on this list. Though it should be the method 
of choice, it has not been used for t he purposes which are t he focus of this paper. An 
example is the modelling of M15 bv iBaumgardt et aD ( 20031 ). who constructed a small 
version which, when scaled up, corresponded approximately to the conditions expected 
for M15. The difficulty i s that unsealed iV-body models are pract ically limited to N of 
order 10 5 at present (e.g. IBaumgardt fc Makmoll2003l: lH~urlevll2007l ). whereas the median 
for the globular clusters at the present time is about 5 x 10 5 . Though it might be hoped 
that the gap would be bridged by the next generation of computers, it must be recognised 
that all globular clusters at the present day have lost substantial numbers of stars, and 
the primordial median must have been higher. Note that it has become possible only 
relatively recently to carry out full simulations of open star clusters. The initial mass of 
M67, for instance, is estimated at about 19 000 Mq, i.e. about 10 times its present mass, 
and thi s simulation, with a realistic complement of primordial binaries, took of order 1 
month (jHurlev et al.ll2005l ). 



1.2.2. The observational constraints 

Whichever method is chosen to model a globular cluster, there are a number of obser- 
vational constraints to be satisfied. In historical order of fi rst use we have 

(a) Surface brightness and/or star counts, starting with Von ZeipeJ (1 19081) : 

(b) Radial velocities, wh ether central averaged values ( Illingworthl 1976 ) or radial ve- 
locities of individual stars (IGunn fc Griffin! 1 19 79) 



(c) Pulsar accelerations ( Phinnev 19931 : Grabhorn 1993) 
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(d) Deep luminosity/mass functions, starting essentially with the advent of studies by 
several authors using HST (1995) 

(e) Accurate proper motions ( van den Bosch et al.l [2006) . (We refer here to models 
built to satisfy constraints imposed by observations of internal proper motions, and not 
to the use of proper motions to establish membership.) 

1.2.3. The Monte Carlo Model 

This paper focuses on an application of the Monte Carlo code developed essentially by 
Giersd (|l998 . 2001 . 2006h . In this approach we assume spherical symmetry and dynamic 
equilibrium, and characterise each star by its energy E and angular momentum J. The 
code repeatedly alters E and J to mimic the effects of gravitational encounters, using 
the theory of relaxation. The same theory underpins Fokker-Planck codes, and the ba- 
sic Monte Carlo code provides essentially a Monte Carlo solution of the Fokker-Planck 
equation. 

The Monte Carlo code is rather suitable for the addition of a number of other process, 
chiefly: 

(a) The galactic tidal field, which is treated as a cutoff 

(&) Binaries, whose interactions are treated using cross s ections ( f rom [spitzerl jl987t ) 
for interactions with single stars, and expressions based on iMikkola (jl983l Il984alfbl ) for 

interactions between binaries) 

(c) Stellar evolution of single stars ( Hurley et al. 2000h and binary stars ( Hurley et al 
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Each of these requires some comment, though further details are given by Giersz & Heggie 
in this volume, and in Sec. 4 below. 

(a) There are significant differences between a tidal field and a ti dal cutoff, as thes e 
lead to somewhat different scalings of the dissolution time with N (jBaumgardtl 1200 if ) . 
We have attempted to mimic this with a mass-dependent lowering of the escape energy. 

(b) The Monte Carlo code described here has no triples, and so hierarchical triples 
(which are a common product of binary-binary encounters) have to be bypassed. Fur- 
thermore, the use of cross sections hinders the inclusion of star-star collisions during 3- 
and 4-body encounters. Finally, the cross sections are not well known for unequal masses . 



(c) Stellar evolution is implemented via the McScatter interface ([Heggie et al.1 12006). 



Besides the stellar evolution packages of Hurley et al (referenced above) it also inter faces 
to the stellar evolution package SeBa in starlab ([Portegies Zwart fc Verbuntlll996l ). At 
present, however, the latter is limited to solar metallicity, which is unsuitable for the 
study of old globular clusters. 

In view of the above approximations and uncertainties, the testing and calibration 
of the Monte Carlo code against the results of TV-body models (in the regime of small 
enough N) is an essential safeguard. Such studies are described by Giersz & Heggie in 
this volume. 



2. The globular cluster M4 

This, one of the very nearest known globular clusters, was selected at the meeting 
MODEST-5 (Hamilton, Canada, 2004) as a target for concerted observational and theo- 
retical effort, but so far little theoretical work has been carried out. Table 1 summarises 
some data for this fascinating object. 

The proximity of M4 makes deep observational study possible. (See, for example, the 
poster by Sommariva et al in this volume.) For theoretical purposes too it is well placed 
for study because its binary population appears to be modest, and its initial mass may 



4 



D.C. Heggie and M. Giersz 



Table 1. Properties of M4 



Distance from sun 


a 1.72kpc 


Distance from GC 


5.9kpc 


Mass 


63 000 M, 


Core radius 


0.53 pc 


Half-light radius 


2.3 pc 


Tidal radius 


21 pc 


Half- mass relaxation time (Rh) 


660 Myr 


Binary fraction 


a l-15% 


[Fe/H] 


-1.2 


Age 


b 12Gyr 


A v 


a 1.33 



References: All da ta are from the current version of the catalogue of IHarrisl (1 19961 ). except a 
Richer ct al. (2004) (though this is not always the original reference for the quoted number) and 
^Hansen et al.l (|2004l ). 



not have been very high, as we shall see. One compli cation for the Monte Carlo code, 
however, is that the orbit appears to be very elliptical (jDinescu et al. I ll999l) . whereas we 
must assume a steady tidal field. 

Table 2 describes the initial conditions which we adopted for this exercise. The primary 
observational data which we attempted to fit were 



(a) The surface brightness profile (ITrager et al.lll 995) 



(b) The radial velocity dispersi on profile (iPeterson et al.l ll995) 

(c) The V- luminosity function ([Richer et al. 20041 . from which we considered the re- 
sults for the innermost and outermost of their four annuli) 

though several other observational comparisons will be described below. We do not have 
a systematic way of arriving at a best choice of initial parameters, though a possible 
approach is described towards the end of this paper. We began with a scaled-up version 
of the models we developed for the old open cluster M67 (see Giersz & Heggie in this 
volume), but found that a binary population of = 50% tended to produce a model 
with too low a concentration. Reducing the binary concentration to 5 or 10% produced 
a satisfactory surface brightness profile, but was somewhat too massive, because of an 
excess of l ow-mass stars, corresponding to a poor fit with the luminosity function. Ac- 
cording to lBaumgardt fc Makind ([2003) it might be possible to correct this by devising 
a model which lost mass at a higher rate, but instead we elected to change the slope of 
the low-mass IMF from the canonical value of a = 1.3 (|Kroupall2007l) to a = 0.9. (There 
is some justification for a lower value for low-metallicity populations.) By some experi- 
mentation we arrived at a model which gave a fair fit to all three kinds of observational 
data; see Table 3, and Figs 1-4. Much of the disagreement in the total luminosity is due 
to our assumed distance to t he cluster, which is significantly smaller than the value of 
2.2kpc given bv IHarrisl (|l996l) . though our surface brightness profile is also a little faint 
on average. The disagreement in the inner luminosity function at faint magnitudes may 
be attributable to the fact that the theoretical result assumes 100% completeness, while 
the observational data are uncorrected for completeness. A typical plot of a completeness 
correction is given bv lHansen et all (|2002L Fig.3). 

It is worth noting that no arbitrary normalisation has been applied in these compar- 
isons between our model and the observations. The surface brightness profile, for example, 
is computed directly from the V-magnitudes of the stars in the Monte Carlo simulation. 

Figure [5] shows the colour-magnitude diagram of the model. This is of interest, not so 
much for comparison with observations, but for the presence of a number of interesting 
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Table 2. Initial parameters for M4 



Fixed parameters 
Structure 
Stellar IMF 

Binary mass distribution 
Binary mass ratio 
Binary semi-major axis 
Binary eccentricity 
Metallicity Z 
Age 



Plummer model 
Kroupa double power law 
Kroupa et al. (1991) 
Uniform 

Uniform in log, 2(Ri + R 2 ) to 50AU 
Thermal, with eigenevolution (Kroupa 1995) 
0.002 
12Gyr 



Free parameters 
Mass 

Tidal radius 
Half-mass radius 
Binary fraction 

Slope of the lower mass function 



M 
rt 
rh 
fb 

a (Kroupa = 1.3) 



Table 3. Monte Carlo and King models for M4 



Quantity 


MC model 
(t = 0) 


MC model 
(t = 12Gyr) 


King model 
(Richer et al. 20041 


Mass (M Q ) 


3.40 x 10 5 


4.61 x 


10 4 




Luminosity (Lq) 


6.1 x 10 6 


2.55 x 


10 4 


6.25 x 10 4 


Binary fraction ft, 


0.07 


0.057 







Low-mass MF slope a 


0.9 


0.03 




0.1 


Mass of white dwarfs (Mq) 





1.81 x 


10 4 


3.25 x 10 4 * 


Mass of neutron stars (Mq) 





3.24 x 


10 3 




Tidal radius r t (pc) 


35.0 


18.0 






Half-mass radius rh (pc) 


0.58 


2.89 







*: this is the quoted mass of "degenerates" 



features. The division of the lower main sequence is simply an artifact of the way binary 
masses were selected (a total mass above 0.2M Q and a component mass above 0.1M Q .) 
Of particular interest are the high numbers of merger remnants on the lower white dwarf 
sequence. There are very few blue stragglers. Partly this is a result of the low binary 
frequency, but it is also important to note that some formation channels are unrepresented 
in our models (in particular, collisions during triple or four-body interactions, though if 
a binary emerges from an interaction with appropriate parameters, it will be treated as 
merged.) These numbers also depend on the assumed initial distribution of semi-major 
axis, which is not yet well constrained by observations in globular clusters. 

Photometric bi naries are vis i ble in Fig [51 and these are compared with observations in 
the inner field of iRicher et al.l (|20041 ) in FigJS) In this figure, the model histogram has 
been normalised to the same total number of stars as the observational one. We made 
no attempt to simulate photometric errors, but the bins around abscissa = -0.75 suggest 
that the binary fractions in the model and the observations are comparable. Fig[7] shows 
that binaries have evolved dynamically as well as through their internal evolution. In 
particular the softest pairs been almost destroyed. 

By 12 Gyr the binaries exhibit segregation towards the centre of the cluster, but 
perhaps in more subtle ways than might be expected (Figs l8l9| ). When all binaries are 
considered, there is little segregation relative to the other objects in the system. (Most 
binaries in our model are of low mass.) But if one restricts attention to bright binaries, 
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Surface brightness 



Velocity dispersion 




1000 10000 



Figure 1. Surface brightness profile of our 
Monte Carlo model , com pared with the 
data of iTrager et~al] (119951 ). 




Figure 2. Velocity dispersion profile of 
our Mo nte Carlo mode l, comp ared with the 
data of lPeterson et all (|1995T ). 



Innermost annulus 



Outermost annulus 





Figure 3. Luminosity function of our 
Monte Carlo model at the median radius 
of the innermost annulus in iRicher et al.l 
(2004), compared with their data. 



Figure 4. Luminosity function of our 
Monte Carlo model at the median radius 
of the outermost annulus in IRicher et al.l 
(2004), compared with their data. 



which we here take to mean those with My < 7 (i.e. brighter than about two magnitudes 
below turnoff), the segregation is very noticeable (Fig[5]), with a half-mass radius smaller 
by almost a factor of 2 than for bright single stars. Still, bright binaries are not nearly 
as mass-segregated as neutron stars (Fig|8]), which, incidentally, receive no natal kicks in 
our model. 

These data do not reveal one very interesting feature of our model, which is that 
it exhibited core collapse at about 8Gyr. Subsequently its core radius is presumably 
sustained by binary burning. Even non-primordial binaries may be playing a role here. 
To the best of our knowledge it has not previous ly been suggested that M4, which is 
classified as a "King" cluster bv lTrager etail (|l993h . is a post-collapse cluster. This raises 
the long-dormant question of how it is possible for some clusters to exhibit collapsed cores 
if they also come with significant populations of primordial binaries. 
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CM Diagram 



10 



15 



20 




-0.5 



Figure 5. The colour-magnitude diagram at 12Gyr. Green: single stars; red: binaries; blue 

pluses: collision or merger remnants. 




Distribution of semi-major axes 



V offset from single star MS 



Figure 6. Histogram of V-oSset from the 
main sequence, compared with the corre- 
sponding data from the innermost annulus 
studied by Richer et al. See text for details. 




100 1000 10000 100000 

Semi-major axis (solar raOii) 



Figure 7. Distribution function of the 
semi-major axes of the binaries at OGyr and 
12Gyr. Units: solar radii. 



3. The search for initial conditions 

Each Monte Carlo model for this cluster takes a few days. Therefore the problem of 
finding appropriate initial conditions is a significant one. One needs a good starting guess, 
and then a rapid method for iterative improvement. Our techniques for dealing with these 
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Radial distribution 




0.1 1 10 

Radius (pc) 

Figure 8. Radial distribution functions at 
12Gyr. Units: parsecs. The key identifies 
the class of object included. The distribu- 
tions of white dwarfs and binaries are al- 
most identical. 



Radial distribution of bright objects 




Figure 9. Radial distribution functions at 
12Gyr, showing the extent of segregation 
between bright single and bright binary 
stars (Mv < 7). Also shown for compari- 
son is the distribution for all objects, as in 
Fig(S] Units: parsecs. 



issues are still primitive, but have evolved in the course of this research. At the iterative 
stage we have employed scaled- up small models (with as few as 10 4 objects sometimes), 
which are designed to relax at the same rate as a full-scale model. This requires a scaling 
of the length-scale, which does violence to binary interactions, but nevertheless it has 
been very useful. In this section we focus on the issue of finding starting values. 

Consider first the problem of forwards evolution. Useful formulae for M(t) are given by 



Lamers et al.l (|2005l h which we have generalised to other IMF's and metallicities. To this 
we have added formulae for rj- l {t) 1 based on very simple notions of adiabatic expansion 
(in response to mass-loss from stellar evolution) and tidal truncation. For the evolution 
of the co re we have fitted simple expr essions to the results on the time of core collapse 
given bv lBaumgardt fc Makinol (|2003lh extended by new iV-body simulations to a wider 
range of initial concentrations. Very simple expressions for r c (t), consistent with the 
initial and final values, can then be employed as a first approximation. Finally we drew 
from Baumgardt & Makino a relation between M and a (the slope of the lower mass 
function.) Putting these relations together, we have constructed a tool which we refer 
to as Quick Cluster Evolution, following S. Portegies Zwart. To apply this to generate 
initial conditions for our simulations, we can run QCE iteratively in reverse. 

Further development of this tool should include the addition of binary heating, which 
certainly influences the evolution of the half-mass radius unless the primordial binary 
fraction is low enough, and more concentrated initial conditions than the King models 
to which we have restricted QCE so far. 



4. Discussion 

It is shown by Giersz & Heggie (this volume) that Monte Carlo models can provide 
similar results to iV-body models, in the range where comparison is possible, with similar 
physics (binaries, stellar evolution, etc.), except for a number of restrictions: 

(a) Use of a tidal cutoff, instead of the tidal field. Though Sec 1 1 . 2 . 3l summarises our 
current approach to this problem, other treatments are possible, and worth trying. 

(b) Use of a static tide. The effects of tidal shocks have been studied by a number 
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of authors (e.g. Kundic k. Ostrikerlll995 ). and it would be possible to add the effects as 
another process altering the energies and angular mom enta of the stars in the simulation. 

(c) Rotation: it has been shown (|Kim et al.l 120041 ) that, to the extent that rotating 
and non-rotating models can be compared, rotation somewhat accelerates the rate of 
core collapse. Rotation is hard to implement in this Monte Carlo model, however. 

(d) Use of cross sections for triple/quad interactions: this limitation could be o vercome 
by direct integration of the interactions, as is done by iFregeau k, RasTol ( 2007 ) in their 
version of the Monte Carlo scheme. It will be important to do so, as it would remove the 
dependence on cross sections which are not well known for unequal masses, and would 
also permit us to inclu de the collisions betwe en stars which commonly occur in long-lived 
few-body interactions ( Hut fc Inagaki 19851 ). 

(e) Neglect o f triples: these are also commonly produced in binary-binary encounters 
( Mikkolalll984al) . and it is desirable to include these as a third species (beyond single and 
binary stars). Their observable effects may be small, but of course there is one intriguing 
example in the very cluster we have focused on here ( Thorsett et al.l[l993l ). 

Despite these limitations, not all of which are easily curable, Monte Carlo models are 
feasible in reasonable time for globular clusters, which are too large for direct iV-body 
models. They yield predictions for mass segregation, luminosity functions, distributions 
of binary parameters, anisotropy, and many other kinds of data, which can hardly be 
obtained in any o ther way. (The only com parable method of which we are aware is 
the hybrid code of Giersz &: Spurzem 20031 .) Even when TV-body simulations eventually 
become possible, Monte Carlo models will remain as a quicker way of exploring the main 
issues, just as King models have continued to dominate the field of star cluster modelling 
even when more advanced methods (e.g. Fokker-Planck models) have become available. 

In addition to some of the possible improvements mentioned above, it is our inten- 
tion to extend the approach to a number of other objects, including a "collapsed-core" 
cluster such as NGC6624. We welcome all suggestions for observational or theoretical 
comparisons, either on M4 or on other objects. 
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